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The paper deals with the numerical solution of the nonlinear 
Ito stochastic differential equations (SDEs) appearing in the unrav- 
elling of quantum master equations. We first develop an exponential 
scheme of weak order 1 for general globally Lipschitz SDEs governed 
by Brownian motions. Then, we proceed to study the numerical in- 
tegration of a class of locally Lipschitz SDEs. More precisely, we 
adapt the exponential scheme obtained in the first part of the work 
to the characteristics of certain finite-dimensional nonlinear stochas- 
tic Schrodinger equations. This yields a numerical method for the 
simulation of the mean value of quantum observables. We address 
the rate of convergence arising in this computation. Finally, an ex- 
periment with a representative quantum master equation illustrates 
the good performance of the new scheme. 

1. Introduction. 

1.1. Objectives. The primary objective of this paper is to develop an ef- 
ficient scheme for the computation of E{Zt,AZt), where A S C'^''^, (•, •) is the 
standard scalar product in C^, and Zt satisfies the fohowing Ito stochastic 
differential equation (SDE) on C^: 
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where ||^o|| = 1; W is an n-dimensional (^5^j)-Brownian motion, 

n 

G = -iH - ^ ^ LlLk, 

■^(^) = ^^^^'\^^^k^)^kZ - ^Re^(2,Lfcx)z), 

fc=i 

provided z G C^, and for any /c = 1, . . . , n and z G 

(3) Ekiz) = LkZ -Re{z,Lkz)z. 

Here L^,- G C'^''^ for all k = 1, . . . ,n, H is a d x d self-adjoint matrix, and 
(r2,5^, P, {dt)t>o) is the underlaying filtered probability space. Our main mo- 
tivation came from the numerical simulation of the evolution of open quan- 
tum systems. To help to shed light on our problem. Section 1.2 below looks 
closely at this application. 

In this work we follow the strategy of constructing exponential schemes 
adapted to the characteristics of (1). This approach is partially motivated 
by the good behavior of the exponential integrators in the solution of certain 
class of stiff ordinary differential equations, for example, those associated to 
both time-dependent Schrodinger equations and oscillatory electric circuits 
(see, e.g., [12, 13] for more details). Another motivation came from a num- 
ber of numerical experiments which illustrate the good performance of the 
exponential schemes for real SDEs with additive noise whose numerical so- 
lution by the standard explicit schemes presents numerical instabilities (see, 
e.g., [4, 16, 24]). 

To gain understanding of the exponential methods, this article starts by 
generalizing the Euler-exponential scheme for SDEs with additive noise pro- 
posed in [24] to the context of SDEs of the form 

(4) Xt = Xo+ f\{s,Xs)ds+ f' a{s,Xs)dWs, 

Jo Jo 

where t G [0,T], Xt takes values in M.'^ and b, a are smooth functions with 
bounded derivatives up to appropriate order. Indeed, adapting the method- 
ology employed in [24], we develop the following numerical method: 

Scheme 1 (Euler-exponential). Let ^o,---,Coi ■ ■ ■ ■ ■ ■ , ^m-1' • • • ' 
(,M-i be independent and identically distributed (i.i.d.) real random vari- 
ables with symmetric law, variance 1 and moments of any order. Then, we 
consider the recursive algorithm 



Vrn + J^iKTm, ^m) " Jb{T„„ Vm)Vm) + 
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where = mT/M, Jb = {dib^)k,j=i,...,d and U = (d, ■ • ■ , C)^- Here Vq is 
independent of for all m = 1, . . . , M — 1. 

In particular, we check that the error between ¥jf{Xx) and 'EifiV^^) has a 
linear behavior when / is a smooth function, that is. Scheme 1 achieves the 
first order of weak convergence. Preliminary numerical experiments suggest 
that Scheme 1 should be useful in situations where the eigenvalues of Jb 
have vastly different sizes and their real parts are nonpositive. In these cir- 
cumstances, both the explicit Euler scheme and the implicit Euler scheme 
present relevant time step restrictions in many cases. It will be interest- 
ing to test more carefully Scheme 1 with theoretical and real-life scientific 
problems. 

Combining Scheme 1 with splitting and projection techniques yields an 
efficient numerical method for (1). To be more precise, we may split the drift 
term of (1) into two components to obtain that there exists a continuous 
semimartingale S.^t^ such that 

(5) Zt = ZT^+ f GZrdr + St,Tm 

for all t G [Tm, Tm+i]. Here Tm = mT/M . Then, solving explicitly, the linear 
SDE (5) leads to 

(6) Zt = exp(G(t - Tm))ZT^ + /* exp(G(t - r)) dSr,T^, 

for any t G \Tm-,Tm+i\- Letting Zm ~ Zt^ and replacing the right-hand side 
of (6) by random vectors with similar first three moment properties, we 

Z A'l 

arrive at the weak approximation ^^^'^ of ^t„+i , where 

(7) = exp(G^) l^z + Diz)I^ + y^E , 

with ^Q, . . . , , . . . , Cm-I' • • • ' ?A/-i Scheme 1. Now, the growth of 

is stabilized by using a projection technique from the numerical treat- 
ment of ordinary differential equations with invariants. To be precise, since 

||Zt|| = 1, we project ^^^^ onto the surface of the unit ball. From this we 
derive the following numerical method: 

Scheme 2 (Version of the Euler-exponential method). Let Zq^ be a 
random variable independent of ^i, . . . ,^m-i satisfying H^q"^!! = 1. Then, we 
define recursively 
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where 




0, 

zl\\z 



if 2; = 0, 
if 



The behavior of Zm is tested by means of a numerical experiment. Some- 
times the apphcation of the projection techniques deteriorates the perfor- 
mance of the numerical method. For instance, Section VII. 2 of [10] presents 
an example where a projection procedure destroys the correct qualitative 
behavior of a symplectic Euler method applied to a deterministic Hamil- 
tonian system. In our case, a numerical experiment with a representative 
quantum system, where the eigenvalues of G have vastly different sizes and 
their real parts are nonpositive, illustrates the very good behavior of Z^- In 
this example, both versions explicit and implicit of the Euler scheme fail. 

A secondary objective of this paper is to study convergence properties 
of the stochastic schemes used for computing 'E{Zt,AZt). In particular, our 
interest is focused on the rate of convergence of 'E{Zm.,AZm). Most of the 
existing convergence theory for numerical methods requires that the co- 
efficients of the SDE be globally Lipschitz. This motivates the increasing 
interest in addressing convergence properties of the numerical schemes for 
more general class of SDEs (see, e.g., [11, 20, 31]). In our case, the SDE 
under, consideration is only locally Lipschitz. In fact, the coefficients of (1) 
have nonlinear grow. To deal with this situation, we modify the standard 
arguments due to Talay [27, 28] and Milshtein [21]. Another difficulty in car- 
rying out our theoretical study is that p has a singularity at 0. To overcome 
it, we take some inspirations in [29] and [9]. 

This paper is organized in six sections. Section 2 is devoted to intro- 
duce notation. Section 3 develops Scheme 1. In Section 4 we construct 
Scheme 2. Section 5 provides the rate of convergence of ¥!,{Zjlj , AZj^) . Sec- 
tion 6 presents a numerical experiment. 

1.2. Motivation. We start with a brief exposition of some basic results 
of quantum mechanics. The states of a quantum system are described by 
elements of an adequate complex Hilbert space (f), (•,•)) and the quantum 
observables are represented by self-adjoint linear operators in i) (see, e.g., [5] 
for more details). In the Heisenberg picture, the evolution of the observable 
A under the Born-Markov approximation is given by the minimal solution 
of the adjoint quantum master equation 



(8) 



— Tt = G*Tt + TtG + Y,LlTtLk, To=A. 



Here Tt, Li 



. . . ,Ln are general linear operators in Pi and 



n 



G = -iH-\Y.LlLk 
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with H self-adjoint operator in P|. The operators Lj, with j = 1, . . . ,n, de- 
scribe the effects of the environment and H represents the Hamiltonian. For 
a fuller mathematical treatment of (8), we refer the reader to [8]. 

The concept of quantum trajectories allows the simulation of the evolution 
of the quantum observables. Let us consider the linear stochastic evolution 
equation on [), 

ft " pt 

(9) Yt = yo+ GYsds + y^ LkY^dB^, 

Jo ^^Jo 

where B is an n-dimensional Brownian motion on the filtered complete prob- 
ability space Qj (5^i)t>o)- Then, in the measurement interpretation of 
the quantum trajectories the stochastic process Zt = lt/||yj|| describes the 
evolution of a system conditioned on continuous observation (see, e.g., [33] 
and the references given there). Furthermore, the mean value of the observ- 
able A at the instant t is given by 'EQ{Yt,AYt). In fact, we may see that 

(10) Bci{Yt,AYt) = {yo,Ttyo), 

under certain assumptions (see, e.g., [3, 15]). 

The stochastic Schrodinger equations (1) allows the description of finite- 
dimensional open quantum systems, for example, g-bits models. Let dim < 
+00. Applying Ito's formula, integration by parts formula and Girsanov's 
theorem, we see that there exists a probability measure P, which is equiva- 
lent to Q, such that Zt satisfies (1) and 

(11) BQ{Yt,AYt) = \\yofBp{Zt,AZt). 
Combining (10) and (11) gives 

{yo,nyo) = \\yo\\^'Ep{Zt,AZt). 

Hence, the numerical solution of (1) leads to the numerical description 
of {yo,Ttyo), which represents the mean value of the observable A at the 
instant t. This procedure has been proposed in the physical literature in or- 
der to overcome the difficulties appearing in the direct numerical integration 
of (8) (see, e.g., [25]). It is worth pointing out that the numerical schemes 
for (8) exhibit serious numerical instabilities and the dimension of the state 
space of (8) grows up very fast to +oo when dim i) — > +oo. On the other 
hand, the computation (yoi'^tyo) by means of the numerical solution of (9) 
presents drawbacks. In this case, a large number of numerical experiments 
show the blow-up of the trajectories of the explicit Euler method even for 
small size of the discretization step. Furthermore, in many situations the 
implicit Euler scheme, defined as in [17, 20, 22, 31], tends to the origin very 
fast. 

Finally, the efficient numerical solution of (1) also plays an important role 
in the study of many infinite-dimensional quantum phenomena, for instance. 
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harmonic oscillators. Let dimt) = +00. Then, we may choose an adequate 
sequence {i)d)d of finite-dimensional Hilbert spaces such that EQ{Yt, AYt) 
is approximated by EiQ{Yt^d, AYt,d) , where Yf^d is the continuous adapted 
stochastic process with values on i)d given by 

Yt,d = PdYo+ / GdYs,dds + J2 / PdLkYs,ddB^, 
Jo ^^^Jo 

with Pd'-t) ^ t)d the orthogonal projection of f) over f)^ and 

n 

Gd = -iPdH - \ J2 PdL*PdLj 
i=i 

(see, e.g., [23], where the rate of convergence of this approximation is stud- 
ied). Thus, similar arguments to those used in the previous paragraph give 
rise to our main problem. 

2. Notation. Throughout this paper we assume that the filtered proba- 
bility spaces satisfy the usual hypotheses (see, e.g., [7, 26]). We will denote 
by Et the conditional expectation with respect to For simplicity, we re- 
strict our attention to the equidistant partitions of the time interval [0,T], 
that is, time discretizations of the form {T^)m=o,...,M, with T^^ = mT/M. 
To shorten notation, sometimes the explicit dependence on the discretiza- 
tion step T/M will be suppressed except where we wish to emphasize its 
role. For example, we will write Tm instead of if no misunderstanding 
is possible. We will use the same symbol K{-) (resp., K and q) for different 
positive increasing functions (resp., positive real numbers) having the com- 
mon property to be independent of M. Furthermore, it is understood that 
q is greater than or equal to 2. 

Let A G C''"^. Then, the symbol will stand for the transpose of A. 
Furthermore, A^'^ will be the {k,j)th component of the matrix A and ||^|| = 

y^SL=i X]j=i l^'^'-' P- For any x,y G C^, we will write {x, y) = J2k=i ^^V^ ^^"^ 

X = {x^, . . . , x'^). 

For each Z € N, we define Vi to be {1, ... , d}'. For any p G Vi, with / G N, 
and xeT'^, with T G {M, C}, we set 

F,^{x) = n x^' 

k=l 

and dPg{x) = ■ ■ ■ -^^g{x), provided that 5 : T*^ — > T is smooth enough. 

The symbol stands for the identity operator and Vq = {0}. We say that 
a family of functions {fe : [0,T] x M*^ ^ ^)9ee belongs to C^([0,r] x M'',M) 
if for any p^^Vi, with I < L, 
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(i) d^f$ is a continuous function whenever 9 £ Q, 

(ii) \dPfeit,x)\ <K(r)(l + for alHG [0,T], xGR'^, and G G. 

Moreover, we write {fe-^'^ ^ ^)eGe G C^(M'^,K) if (/e)0Ge satisfies the con- 
ditions (i) and (ii). 

3. Euler-exponential scheme for general SDEs. To shed some new hght 
on exponential schemes, this section develops an exponential method for (4). 
First, we introduce a local exponential representation of X. Second, we de- 
rive a first weak order exponential scheme. Finally, we deal with the conver- 
gence analysis of the new scheme. 

3.1. Euler-exponential scheme. Throughout this section, we assume that 
b and a satisfy the standard conditions for the existence and uniqueness of 
strong solutions of (4) (see, e.g., [1, 26]). Furthermore, for any s € [0,r] and 
j; E R'^, we consider the adapted stochastic process X^'^ defined by 



(12) Xp^ = x+ f bir, X,^'") dr + /* a{r, X, 

J S J S 

for all t G [5,r]. For abbreviation, we set :=Xi"'^. 

The following theorem states a local representation of X of exponential 
type. 

Lemma 1. Let d^b be a continuous function for each p G "Ps. Suppose 
that dtb and dtd^b, with k = 1, . . . , d, are also continuous functions. If ^ is 
a ^T^n-i^CLndom variable taking values in W^, then for any t G [Tm,T„-i^i], 

JT-rn 

(13) + /* e''^^-'^^^'-'\b{T„,,0 - Jb{T„,,00ds 

7' L{b){u,Xl-^^)dv)ds + Rt,T^. 
Recall from Section 1 that Jb{Tm,^) = {dlb^{Tm,0)k,j=i,...,d- addition, 



^ T ' 

'If 



* Jh{Tm,i){t-s) 



and Rt,Tm is given by 

f e-^''(^-'«)(*-^) r (r Li{dib){u,X^"^'^)dWl\V{r,X^^'^)drds 



j=l i=Q • 
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d n f 



+EE 

j=l i=0 



r (f U[dih){u,Xl-'^)dWl)a^^\r,Xj-'^)dWrds, 



where 



k=l k,l=l 



and for any i = I, . . . ,n, 

k=l 

Proof. Using Ito's formula (see, e.g., [6]), we obtain, for any s > T^, 

b{s,X^^'^) = b{Tm,0+ f L{b){u,X^^^^)du 

JTrn 

+ X: r dib{u,x^-'^)d{x^-'^y. 

Substituting this result into (12), we see that, for any t > Tm^ 
X^"^'^=^+ f b{T^,Ods+ f (f L{b){u,X^-^^)du\ds 

+ E r (r dib{u,X^-'^)d{X^-'^y)ds+ f a{s,X^-'^)dWs. 
Then, applying Ito's formula to each 5^6, with z = 1, . . . , d, yields 
(14) = ^ + /* JHTm,OXl""^ ds + St,T^ yt G [Tm, Tm+l] , 

where 

St,T^= f {b{Tm,i)-Jb{Tra,mds+ f a{s , X^^^^) dW s 
+ f (f Lib){u,X^-'^)du\ds 

JTfrt \JTrn J 

+ f if if U{dih){u,Xl-^^)dWl)V{s,X^-^^)ds)dr 

+ / / Hdib){u,X^-'^)dW:)a^'-{s,X]'"^'^)dWs)dr, 

^ — 1 ^ — r\JTm \JTm \JTm / / 



i=ii=o- 
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with := u. Since S.^Tm is a continuous semimartingale, the hnear SDE (14) 
has the exphcit solution estabhshed in (13) (see, e.g., [26]). □ 

We are now in position to deduce an exponential scheme of weak order 1 
for (4). Since -^t™+i =^t™+^^'"' start by replacing ^ by Xt^ in (13). 
Then, we neglect the terms of the last line of (13) because they involve 
multiple integrals. Applying the Ito-Taylor formula to a{s,Xs), we see that 
a{s,Xs) can be approximated by a{Tm,XT^) in the first line of (13). Fur- 
thermore, we substitute Xt^ for Xt^ in (13), where Xt^ is so chosen that 
it approximates Xt^ in a weak linear sense. From this we may conclude 
that, for any t £ [Tm,Tm+i], 

Xt^yf^"^'^, 

where for each J^^^^ -random variable ^, 

y/'*^ = exp(j6(r„,,e)(t-r™,))c 

(15) + r exp{Jb{T^,0{t-s)){b{T^,0-JHTm,OOds 

+ f exp{Jb{T^,^){t-s))a{Tm,C)dWs. 

X A'l 

In order to replace Yrp'^"'^' by other random variables with similar first 
three moments, we use arguments similar to those in Section 3.2 of [24]. In 
fact, we approximate the integral 

rTm+i 

exp( JKT^, Xt,„ ) {t - s)){b{T^, Xt^ ) - Jb{T^, XtJXt„^ ) ds 

by a classical rectangle. Finally, we look for a linear function HMiTmrr) 
such that 

with 1^0) •• • )^M-i defined as in Scheme 1, is the approximation of 

given by a classical rectangle rule. This yields Scheme 1 defined in the In- 
troduction, that is, the method 



Vm+i = eyi-g(^Jb{Tm,Vm)j^ 



^[Vm + ^{b{Tm,Vm) " Jb{Tm,V^)V^) + ]l ^^(^m, Vn,)Cn 
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Remark 1. Let the eigenvalues of the Jb have nonpositive real part. 
Since 

exp(j6(r^,y„)^) = (^i-jb{T„M^^ \o(^(^Jb{Tm,v^)^y^, 

Scheme 1 leads to the following version of the implicit Euler scheme which 
avoid the solution of nonlinear equations systems: 

Remark 2. To carry out the computation of exp( J6(Tm, Vm)T/M) times 
a vector u in the implementation of V, we may use Krylov approximations 
with Lanczos process (see, e.g., [12, 14]). In fact, many numerical exper- 
iments with high-dimensional problems illustrate the good performance of 
this numerical method. Furthermore, Hochbruck and Lubich [12] proved that 
the convergence of Krylov methods for exp( J6(Tm, Vm)T / M)u is faster than 
that for the solution of the linear equation (I — Jb{Tm, Vm)T/M)x = u, which 
is required in both methods and the usual implicit Euler scheme when 
b is linear. Alternative methods are Fade approximations, Strang splitting, 
Chebyshev approximations and Magnus integrators. 

3.2. Rate of convergence. Similar to the Euler scheme (see, e.g., [2, 32]), 
the error between E/(Xt) and E/(V^^) can be expanded in powers of T/M 
under general enough conditions. In particular, adapting the arguments used 
in [24] for studying the rate of weak convergence of the Euler-exponential 
scheme for SDEs with additive noise, we obtain the next theorem. 

Theorem 1. Let E||Xo||'' < +oo for any q > 2. Assume that b and 
a are continuous functions such that d^b and d^a are bounded continuous 
functions for every p&Vi, with / = 1 , . . . , 9 . Furthermore, suppose that the 
components of f^d^b, ^Oga and dl§-^a belong to C°([0,r] x ]R^,M) 

whenever p£ Vi, with I = 0,1,2. Let the components of -^pb, belong 
to C°([0,T] X M'^,M). // / eC]^(M'^,M), then there exists a continuous func- 
tion with (^'(s,-))se[o,T] GCp(]R'^,M), such that 



(16) 



Ef{XT)-Ef{Vif) -Lj\^,^s,Xs)ds 



<K(T)(1 + E||X, 



ol 



r\2 



provided that Vq have the same distribution as Xq. 
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Before we prove Theorem 1, we present a series of observations and re- 
sults. We start by considering the probabihty space (j7,(J5,P) that arises 
from the completion of the product measure space induced by (0,3^,P), Vq, 
and the random variables ^q, . . . , ^a/-i- Then, we combine the completion of 
{"St c(Vo^^, ^fc : k < [tM/T] — l))t>o with a limit procedure to construct the 
filtration {&t)t>o that satisfies the usual hypotheses (see, e.g., IV 48 of [7]). 
Let us use from now on the same letter to designate a random variable and 
its natural extension to the Cartesian product space Cl, for instance, {Wt)t>o 
also denotes the stochastic process {Wt o PrQ)t,>o, where Pr^ is the projec- 
tion of ft over Q. Thus, W is an n-dimensional ((J5j)-Brownian motion, Xq 
and Vq are 6o-™easurable, and for any m = 0, . . . , M — 1, is both 0t„+i- 
measurable and independent of (St^. Therefore, we only need to verify (16) 
for X and V defined on {n,&,P, {(St)t>o). 

Lemma 2 recalls well-known results. They may be deduced using Ito's 
formula, the existence of a smooth version of the stochastic flow x i— > X^'^ , 
and induction (some details may be found, e.g., in [18, 30]). 

Lemma 2. Fix /3 S NU {0}. Suppose that b and a are continuous func- 
tions such that d^b and d^a are bounded continuous functions on [0,T] x 
for allpeVi, withl = l,...,P + 2. Let {ge)e&e belong to C^+2(M'^,M). Set 

ue{s,x) = ¥.{ge{XT)/Xs = x) = E(ge(^r'")), 
whenever s G [0,r]. Then {ue)eee e C^+^(R'^,M) and for all 6* G G, 
d 

—ug{s,x) = —C{u0){s,x) ifs£[0,T] andx^R'^, 

(17) 

ug{T,x)=ge{x) i/xGM^ 

where 

C = j2b'dl + \j:{aa^rdt^K 

k=l k,l=l 

Furthermore, for each p&Vi, with I = 0, . . . , (3, we have ^^fug is a contin- 
uous function and 

d ^ 

—gp^UB = -dlC{ue). 

In the sequel, the symbol also denotes the conditional expectation with 
respect to 0t. 

Lemma 3. Let t£ [Tm,Tjn+i]- Then, for any (5t^ -measurable random 
variable ^, we have 

(18) BT„Mt,yt^)=9{Tm,0 + ET„J^ (^^g(s,y/) + £T„,c(5)(s,i;^))ds, 



12 C. M. MORA 

with 

d 



Cr,d9){s, x) = Y^iJbir, Ox + 6(r, - Jb{r, d^Ms, 

k=l 



k,l=l 

provided that all of the derivates of g:[0,T] xW^ ^M. appearing in (18) are 
continuous. 

Proof. Observe that 

JTfji 

+ r a{Tm.,OdWs. 

Then we use the Ito formula to obtain (18). □ 

The proof of Lemma 4 is based on the Burkholder-Davis-Gundy inequal- 
ities and the discrete Gronwall-Bellman lemma. We omit it because it may 
be proved in much the same way as Lemma 4.3 of [24]. 

Lemma 4. Let assumptions of Lemma 1 hold. Then 

(19) e( sup ||yi^r)<i^(r)(i + E(||yo"'r)) 

Vm=0,...,M / 

and 

(20) E^JIIF^^, - e'K^-y^'^^/'-^Vt'r) < K{T) { + \\Vi'r). 



M J 

In addition, for any '&Tm''^'^''^dom variable taking values in M'^, we have 



etJ sup r/'*^-e^^(^-'«)(*-^™)cr) <K(T)f^y^V+iiein 



(21) 



M 



The next result provides a uniform bound of weak order 1 for the weak 
error between and Vm, with m = 0, . . . , M . 

Proposition 1. Suppose that b and a are continuous functions such 
that d^b and d^a are bounded continuous functions on [0, T] x M'^ for any 
p(zVi, with I = 1, ... ,4:. Assume that the components of db/dt and da/dt 
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belong to C°([0,r] x M'^,M). Let {ge)e&e G C^(M'^,M). // Vq have the same 
distribution as Xq, then 

lEgeiXrJ - Ege{V^')\ < K{T){1 + E||Xor)^ 
for alio £@ and m = 0,...,M. 

Proof. We first decompose the global error as a sum of terms involving 
the solution of a parabolic partial differential equation (17). This method- 
ology goes back to Talay [27, 28] and Milshtein [21]. More precisely, let 
ug{s,x) = Egg{X^^). Then 

m 

iBgeiXrJ - BgeiV^)\ < J] \B{Hi)\ + |E(i?|)|, 

i=i 

where = uo{Tj,Y^^^^) — ne(Tj_i, Ij-i) and 

Hi =ue{T,,V,) -ue{Tj,Vj^,) + ue{T,,V,^i) -ue{T,,Y^j-'), 

with Vj_i = e^p{Jb{Tj_i,Vj-i)T/M)Vj-i. 

Due to Lemma 2, we can use Taylor's formula to obtain 

1=1 ' peVi 

+ RjiV,)+R,{Yf;~'), 

with 

where Hp- j ti x d is a diagonal matrix whose components belong to [0, 1] . As 
in the proof of Theorem 4.1 of [24], using the Cauchy-Schwarz inequality 
and Lemmas 2 and 4, we now see that 

(22) mi\ < K(T)(E||y,_i||^ + 1) (^^y. 

Due to Lemma 2, applying Lemma 3 gives 

^T,,,H'2 = / ' ET,_,(-£(n,)(s,y/^-^) + £T,_i,y,_i(n,)(s,y/^-^))ds. 

JTj_i 



14 C. M. MORA 

Then, combining again Lemmas 2 and 3 yields 



(23) 



Ci{ue){s,YY'''))dsdt 



r f BT^_,CT,_„v,_A^ue))is,Ys^'-')dsdt, 



where 



Hence, (21) and Lemma 2 lead to 

(24) \EHi\ < K(r)(E||y,„ir + 1) (^)'. 

From (19), (22) and (24) we deduce the assertion of this proposition. □ 

We are now in position to show the main result of this section. 

Proof of Theorem 1. Let u{s,x) = 'Ef{X^^). As in the proof of 
Proposition 1, we have 

M 

^fiVu) = E ^ (i?r + ^2") + En(0, Fo), 

m=l 

where Hf = tx(T„, Y"^-^) - u{Tm-i,Vm-i) and 

= u{Tm, Vra) " u{Tm., Vm-l) + u{Tm., Kn_l) " n(r^, Y^"^), 

with Vm-l = exp{Jb(Tm-i,Vm~i)T/M)Vm-~i. 
It follows from Taylor's formula that 

1=1 ''■ p&Vi 

+ Rl{Vm)+Rm{Y^Z-'), 

with 

Rm{x) = ^ E d^u{Tm,Vm-l + Ej^^^{x){x - Vm~i))Fj;{x - Ki-l), 
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where Hp-m d x d is a diagonal matrix whose components belong to [0,1]. 
We may check that, for any p E Pi, with I = 1, ... ,5, there exist a constant 
q and functions Cp- in Cp{W^, M) such that, for all m = 1, . . . , M, 



Vm-1,M 



T 

Kn-l) -Cp-(K^_l)( — 



<K{T){l + \\Vm-i\n 



m) ■ 



Proceeding similarly to the proof of (22), we obtain 



n 



cJVr 



m—1 1 



<KiT){nVm-i\\' + l){^^y 



Thus, the mean value theorem leads to 



(25) 



2 5 



1=1 pePi 



dPu{Tm,V„i-i] 



<i^(T)(l + E||K;„_ir) 



n 



'Cp{Vr 



m—1 J 



As in the estimation of in the proof of Proposition 1, Lemmas 2 and 3 
imply that (23) holds with u instead of ug. Then, due to Lemma 2, applying 
Lemma 3 yields 



ETjm 



T 



2M 



m — 1 — 1 Tjji — I 



^m{r) dr ds dt, 



where A{s,0 = - Ci{u) + C^{u) - 2Cs,d^{u))){s,0 and 

^.m{s) = (2£i(£(n)) +£(£i(n)) - C^'iu) - £2(n))(s, r/-^) 

+ (3£t_„z_,(/:'(^)) - 3£T_,y_,(A(u)))(s,y>-: 



Here 



^-i:(5^')«j+^E(S(-^)")*'- 



It follows from Lemma 2 that 

rp \ 2 



(26) 



E( 



M 



<K(r)(E||K 



m— 1 1 



+ 1^ 



ji \ 3 

M, 
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Let 

5 

^{s,x)=A{s,x)/2 + Y^ ^ dPu{s,x)cj^{x)/V.. 

1=1 p£Vi 

Using Lemma 2, we get that (^(s, •))se[o,T] ^ Cp(M'^,M). Then, it follows 
from Proposition 1 that 

(27) |E(*(r„_i, y„_i) - M/(r„,_i,XT_ J)| < K(r)(E||Xor + 1)^- 



Hence, combining Ito's formula, (25), (26) and (27) give 



<K(r)(E||Xor + i) 



T Y 

M ' 



which completes the proof. □ 

Remark 3. According to Theorem 1, we have 

\Ef{XT) - 2Bf{V,fl) + Bf{Vlf)\ < K{T){1 + E\\Xo\n ' 

provided the hypotheses of Theorem 1 hold. This yields a second weak order 
scheme based on the extrapolation Scheme 1 (see, e.g., [32]). 

4. Euler-exponential scheme for stochastic Schrodinger equations. We 

now turn to our main problem. To be more precise, this section provides an 
heuristic deduction of a version of the Euler-exponential scheme adapted to 
the characteristics of (1). 

The following lemma discusses the existence and uniqueness of the solu- 
tions of (1). 

Lemma 5. Let s > 0. Suppose that ^ is a '^g-random variable with E||^|p < oo. 
Then there exists a unique global continuous solution of the SDE 

ft " ft 

(28) Zt'^ = ^+ / {GZ:^^ + D{Z:.'^))dr + Y^ / i?fc(Z,^'«) dVF,^ 

for allt> s. Moreover, \\Z^'^\ = 1 a.s., provided that \\^\\ = 1 a.s. Recall that 
D is given by (2) and Ei, . . . ,En are defined by (3). 

Proof. Since the drift coefficient of (28) and E^, with k = 1, . . . ,n, are 
locally Lipschitz, applying the truncation method, we obtain that (28) has a 
unique local solution (see, e.g., [19, 26]). That is, there exists a stopping time 
Q such that (28) has a unique solution up to Q. This solution has continuous 
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paths a.s. and limsup^^^^^. ||.^f'*'^|| = oo a.s. on {(^ < oo}. By Ito's formula, 
we have 

" /"TjvAi 

(29) ||Z^i^,f = ||ef + 2^ / Re{Zp^,L,Z:.'^){l-\\Z^^')dW,\ 

where r^v is the first exit time of Z^'^ of {x: \\x\\ < N}. This yields 

E||z;;«^,f ^Elief. 

It follows that Q = +00 a.s. (null set depending on ,^). Hence, there exists a 
unique global continuous solution of (28). 

Let St = Efc=i /i Re(Z^^'«, LkZp^) dW^!^. Then iSt)t>s is a continuous semi- 
martingale and 

(30) wz^'^f = + j\i - wzp^f) dSr. 

Thus, the last assertion of the lemma follows from the uniqueness of the 
solution of (30). □ 

We split the drift coefficient of (1) into GZg and D{Zs). Then, analysis 
similar to that in the proof of Lemma 1 shows that, for all t G [Tm,Tm+i], 

ft 



^r,„,z^„ = exp(G(t - Tm))ZT^ + / exp(G(t - s))D{Z 
31 

" ft 

+ / exp(G(t - s))ii;,.(Z^, Zt„J dT^,^ 



Let Zm be a linear weak approximation of Zt^ satisfying ||.^m|| = 1- We now 
approximate Zs'"' ^"^ by Z^ in the right-hand side of (31) to obtain 

Zi«e«(*-^-)z„+ r e^(*-^)l)(Z„)fis + X: f e""^'"'^ Ek{Zm) dW^ 

(32) 

for ah t e [r,„,rm+i]. 

Since our goal is to compute 'E{Zt,AZt) , we should approximate the mea- 
sure induced by the right-hand side of (32). To this end, we use the procedure 
employed in Section 3.1 to yield 

where <I> is given by (7). Finally, to include the information that \\Zt\\ = 
1, Zm+i is projected onto the manifold {z S C^: ||2:|| = 1}. We thus get 
Scheme 2 defined in the Introduction, that is, the method 



7^ -T, 
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where 

0, ifz = 0, 



p{z) 



z/\\z\\, ifz/0. 



Remark 4. As in Remark 1, Scheme 2 leads to the following version of 
the implicit Euler scheme. 

Scheme 3. Let /q^ be a random variable with ||/o^|| = 1. Then we set 

5. Rate of convergence. In this section we focus our interest on the 
proof of the following theorem which establishes the linear convergence of 

Theorem 2. Suppose that, for all B £ C^'^, 

(33) \E{Zo,BZo)-E{Z^',BZ^')\ < \\B\\KiT)^. 
If the law of has compact support, then 

(34) \B{ZT,AZT)-B{Zii,AZ^)\<KiT)^. 

The proof of Theorem 2 starts with bounds for the concentration of As 
in [9, 29], the assumption that has compact support allows us to obtain 
this kind of estimate. 



Lemma 6. Let have compact support. Then there exists an increasing 
positive function K2 such that \\ < K2{T), whenever m = 0, . . . , M — 1 

and \\z\\ = 1. Furthermore, there exist 5 g]0, 1[ and a strictly positive con- 
stant Ki independent of both T and T/M such that ||$^_|_]^|| > Ki for all 
T/M < (5, m = 0,...,M- 1, and zgC^ with \\z\\ = 1. 



Proof. Without loss of generality, we can assume that the support of 
.^0 belongs to the interval [—a, a]. Hence, for any z S C^, with ||2;|| = 1, we 
get 



(35) 



M 



D{z) + 



Ft " 

' k=l 



3T 



rp n 



< V ||Lfc|P + 2aW — y IIL 



fc=i 



M 



k=\ 
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For any z G C*, Re(Gz,z) < 0, and so (exp(Gt))t>o is a contraction semi- 
group on C^. This yields 



I m+l II — 



fc=i 



<i^2(r), 



for all z^'Cr with ||z|| = 1. 

Let (^(z) = ||exp(G)z||. From (35), we see that there exists b G]0, 1[ such 
that 



for TjM < 5 and z £ 



+ D{z)- + \l-Y.E,{z)i 

k=l 

1. Thus, 



M 

Z'^ with lUI 



1 

^2' 



(t>[z + D{z)— + \l—Y,Ek{z)ii]>K^, 



k=l 



provided that ||z|| = 1 and T/M < 6. This gives the last assertion of the 
lemma since (p{x) < || exp(Gr/M)x|| whenever T/M < 1. □ 

Similarly to the proof of Theorem 1, in the sequel we consider the com- 
plete probability space (0,6, P) induced by the random variables Z^^ , 
^0)---)^M-i- In addition, (Jl, C5, P, (6f)t>o) will be the filtered probability 
space satisfying the usual hypotheses induced by {Q,<5,P) and the filtra- 
tion {a{ZQ^ '■ k < [tM/T] — l))t>o. By abuse of notation, we use the same 
symbol for the conditional expectation with respect to both and 0t. 

The role of the local approximation Y in the proof of Theorem 1 is played 
here by ^ given by 



ft " ft 

^t'' = z+ / {G^P' + D{z))dr + J2 / Ek{z)dS, 



k 
r ) 



where t>s and for any k = 1, . . . ,n, 



M-l 



[{r). 



m=0 



The next two lemmas provide information about the behavior of ^ . 

Lemma 7. Let have compact support. Then there exists an increas- 

rpM ^ 

ing positive function K4 such that H^f™' || < K^iT), whenever \\z\\ = 1, 
m = 0, . . . ,M — 1, and t G [TJ^ ,T^_^i]. Moreover, there exist A g]0, 1[ and 
a strictly positive constant independent of both T and T/M such that 
-.r«,.„ ^ r^>..r . . . 1^ [T^^T^^J, and 



zeC^ with \\z\ 



1. 
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Proof. For any t G [Tm,Tm+i], we have 

' k=l / 

Hence, analysis similar to that in the proof of Lemma 6 shows the assertion 
of the lemma. □ 

Lemma 8. Suppose that K^i^Kj^ and A are as in Lemma 7. Assume that 
/gC2'4'4([o,T] X 5 X S,C), with S = {z £ : KsiT) < \\z\\ < i^4(T)}. Let 
have compact support. Then for all t G [Tm,Tm+i[, 



(36) =f{Tm,z,z) 

provided that T /M < A and z with \\z\\ = 1. Here 
£i(/)(s,x,y) = ^(^^(s,x,y)(Gx + D(2))'= + ^(s,x,y)(Gy + :DM) 



k=l 

Furthermore, 



= f{Tm,Z,z) 

Trr. 



(37) + ^ (^|^(s, M/j™'^ ) + £.(/)(«, ^r-^ )) 



where ||0/(z,r/M)|| <K{T){T/Mf and 

n d 



1 (7^ f 

£,(/)(5,x,y)=4(/)(s,^,y) + :^E E ^M^,y)E,{zfE,{zy 



+EE(a^(-'-'y)^.(-)'^.(-) 
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Proof. Applying the Ito formula for a general semimartingale, we ob- 
tain, for any t £ [Tm,Tm+i], 

+ /{T„+,}(^)ET„Y/(^,^^•^^)-/(^,^^^'^^) 



- J^^d^^*^ "^^-'^^ M/fr'^ - ^fr'^)" 

- E ^(*' '^^?^)(^^^- ^?^)') • 

This gives (36). Furthermore, expanding 

in powers of (^?-;^^ - ^tZ'^-)' ^nd (^?;;:;^^ - ^?:;%)^ with j = l,...,d, 
we obtain (37). To this end, we combine (36), Taylor's formula and the mean 
value theorem. □ 

Note that the coefficients of (1) are not globally Lipschitz. To over- 
come this difficulty, instead of using the solution of the usual (in this con- 
text) partial differential equation associated to (28), we employ the function 
v:[0,T]xV-^C described by P = {(x, y):x,y£ C^, {y, x) ^ 0} and 

v{s,x,y) = {y,TT-sx)/{y,x), 

where Tt is the solution of the backward quantum master equation (8) with 

Proof of Theorem 2. Let a:[0,T] x x be given by 

a{t,x,y) = {y,TT^tx). 

According to (8) and Ito's formula, we have Ea(t, Zp^, Z^'^ ) = a(s, z, z), for 
any t £ [0,T]. Hence, ^{Z^^ , AZ^^) = {z,tt-sz)- We thus get 

v{s, z, z) = E(Z^'^, AZ^^), 

provided that \\z\\ = 1. Therefore, 

(38) B{Zij,AZ§)-B{ZT,AZT) = Bv{T,Z^,'z^)-Ev{0,Zo,Z^). 
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Suppose that T /M < min{5, A}, with 6 and A as in Lemma 6 and Lemma 
7, respectively. Since ||Zo|| = 1, from (8), (33) and (38), we conclude that 

M 



m=l 



where 



and 



m—1 



We proceed to estimate From the construction of Scheme 2 and a 
simple computation, we see that, for any z G satisfying ||z|| = 1 and p&Vi 
with I = 1,2, 3, we have 

\^T^-.F,{^'rA'' - e^^"-'z) - ^T^_^F,{^IZ-' - e^^l^z)\ < K{T) {^J . 
Hence, combining Lemma 6 with the deterministic Taylor formula gives 

(39) 



|EFn<i^(T)(^^j . 

This follows by the same method as in the estimation of ff™ in the proof of 
Proposition 1. 

It remains to estimate -^2"- According to Lemma 8, we have 



(40) 



Tm-1 



+ Et_, / £^_^(t;)(s,M/^'^-\M/^'^'"-^ ) ds 



+ Zm-l; 



T 
M 

We may now apply Lemma 8 to the terms of the right-hand side of (40) to 
obtain 

'r\2 



(41) 



I'EH^l <K{T) 



MJ 



To be more precise, fortunately a very long computation shows that 

n 

Cz{v){s,Z,z) = {z,TT-sGz) + {Gz,TT-sZ) + ^{LkZ,TT-sLkz) , 

k=l 
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whenever ||z|| = 1. Since (8) leads to 



(42) 

1 



{y,x) 



{y, G*TT-sx) + (y, TT-sGx) + ^(y, L*TT-sLjx) 



we deduce that, for ||z|| = 1, 

dv 

— {s, z, z) + L;,{v){s, z, z) = 0. 
Therefore, Lemma 8 yields that '£1x^-1^2^ equal to 



0„ I) + E.„_, |^(., C-^-- . ) iris 



Trr 



dv 



Hence, (41) follows from (42) and Lemma 7. 

We conclude from (39) and (41) that (34) holds for T/M <min{5,A}. 
Hence, our claim follows from \\Zt\\ = 1 a.s. and ||-^^/|| = 1- O 

Remark 5. We expect that an expansion similar to (16) holds for the 
error Ei{Zt,AZt) - B{Z§,AZl;}). Nevertheless, the proof of this result is 
still in progress. 

Remark 6. We now turn to (9) with dimf) = +00. It is relevant to 
characterize the global error \'Eiq{Yt , AYt) — ^{Zj^j , AZj^j)\ in function of 
M and the dimension of f)m- A step toward this goal was given in [23], where 
the rate of convergence of B{E^^,AE^^J to Bq{Yt,AYt) is studied. Here 
E denotes the numerical solution of (9) by the Euler scheme. An objective 
of this paper is to advance toward the solution of this problem. 

6. Numerical experiment. This section illustrates the performance of 
the scheme Z. To this end, we consider the following representative exam- 
ple of forced and damped quantum harmonic oscillator in the interaction 
representation. 
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Example 1. Returning to Section 1.2, we choose t) = Let {ipk)keZ- 

be the canonical orthonormal basis on the space Then, the domain 

of the operators and a is {j; G ^^(^+) ^ Sfc>o ^l^fcP < +00} and for ah 
m G Z+, a^ifm = \/m + It/Jm+i and 



The Number operator is defined by = a'^a. 

We now simulate the Hamiltonian as H = i{a^ — a) + N. Furthermore, we 
set Li = 0.2a, L2 = O.Ola^, L3 = O.IA^ and L4 = O.la^ 

In Example 1, describes the state space of a single mode of a quan- 
tized electromagnetic field. The operator , respectively a, is the creation 
operator, respectively annihilation operator. Then, for instance, the term 
i{a^ — a) describes a linear pumping and Li simulates the damping due to 
photon emission. 

To test the scheme Z, we set T = 100 and Yq = (pQ. Moreover, we choose fi^ 
as the linear manifold spanned by {ipj : < j < d} with d = 50. The objective 
is to describe numerically EQ(Yt^5o, A^yi,5o) for t G [0,T]. As we comment in 
Section 1.2, this task leads to solve (1) for d = 50. The parameters selected 
allow us to obtain the "true" value of EQ(Yj^5o, A^Vj,5o) by means of the 
solution of the backward quantum master equation (8) associated to our 
SDE. To this end, we use its explicit solution. It is worth pointing out 
that the numerical solution of finite-dimensional backward quantum master 
equations presents serious drawbacks when the dimension of the state space 
is high (see, e.g., [25]). In fact, some of these problems can be observed in 
our example in case d = 100. 

In the numerical experiment, we compare Scheme 2, Scheme 3 and the 
following version of the explicit Euler scheme: 



'^'^ I -E'fc+i/ll-^fc+ill, if -E'/t+i / 0, 
where E^+i = Ek + {GEk + D{Ek))T/M + yT/M E"=i Ej {Ek)Cl ^ith , • • • , 



^0 ; • ■ • )^Af-i' • • • '€a/-i ™ Scheme 2. In all codes, assumes values ±1, 
each with probability 1/2. 

Figure 1 shows the "true" solution and the approximations obtained by 
the numerical schemes. Moreover, Table 1 looks at the dependence of the 
errors eo to the time step size T/M, where 




if m = 0, 
if m > 0. 





2-10-* 



ej(x,M) 



max^^E(Z„iVZ,)-2-10-4 ^ ixfli/iooi^k), Nxfii/.^oi^k)) 




Q I I I I I I I I I I I 

10 20 30 40 50 60 70 80 90 1 00 

1 

Fig. 1. Dotted line: "true" solution, solid line: (a) explicit Euler scheme, (b) implicit 
Euler scheme and (c) Euler- exponential scheme. 



whenever x denotes the numerical method and J G {0,...,100}. Indeed, 
Table 1 presents estimated values of Eq and A. Here A is the maximum of 
the length of the 90 percent confidence intervals taken over the instant of 
times {0, . . . , 100}. We use the batch method to estimate these intervals (see, 
e.g., [17]). 

Table 1 

Errors versus step sizes for the explicit Euler method E, the 
version I of the implicit Euler method and the 
Euler- exponential method Z 



M 


2000 


4000 


8000 


16000 


e{E,M) 


46.6545 


46.7107 


46.6381 


23.5562 


A{E,M)/2 


0.023207 


0.13302 


0.31203 


0.28929 


s{i,M) 


6.6179 


5.5739 


3.9754 


2.5181 


A(/,Af)/2 


0.030248 


0.045375 


0.037721 


0.054798 


e{Z,M) 


0.33533 


0.2236 


0.11426 


0.037446 


A(i,M)/2 


0.066711 


0.059289 


0.077666 


0.098786 




It can be seen from both Figure 1 and Table 1 that the Euler-exponential 
scheme presents a superior performance than the other two numerical meth- 
ods for this example. For instance, the error induced by Scheme 3 with 
M = 16 • 10^ is substantially greater than the error induced by Scheme 2 
with M = 2 • 10^. Furthermore, the accuracy of Scheme 2 is very good for 
large time step sizes. This suggests that Z shows great promise for the long 
time integration of stochastic Schrodinger equations. 

Finally, Figure 2 shows precision-step size diagrams. In particular, this 
figure gives the errors ej(Z,M), with J = 0,50, 100, versus the step size 
T/M . Moreover, it presents the best least square linear approximation of 
each ej{Z,-). From Figure 2, we see that the errors induced by Z closely 
follow a straight line. This is in a good agreement with Theorem 2. 
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